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We demonstrate that frustrated exchange interactions can produce exotic 3D crystals of vortex strings near 
the saturation field (H = Hsat) of body- and face-centered cubic Mott insulators. The combination of cubic 
symmetry and frustration leads to a magnon spectrum of the fully polarized spin state (H > Hsm) with degenerate 
minima at multiple noncoplanar Q vectors. This spectrum becomes gapless at the quantum critical point H = 
//sat and the magnetic ordering below //sat can be formally described as a condensate of a dilute gas of bosons. 
By expanding in the lattice gas parameter, we find that different vortex crystals span sizable regions of the phase 
diagrams for isotropic exchange and are further stabilized by symmetric exchange anisotropy. 


PACS numbers: 75.10.-b, 75.10.Jm, 75.45.+j, 75.70.Kw 

Topological spin structures are of great potential interest 
in future applications of spin-electronic techniques ||T]. The 
skyrmion crystals discovered in the /120-structure metallic 
alloys MnSi and Fei_j^Co^ IE}® and in the Mott insulator 
CuzOSeOs El 0 are prominent examples. While the emer¬ 
gence of crystals of topological structures is reminiscent of 
the Abrikosov vortex lattice of type-II superconductors |17]|8l, 
their origin is completely different in magnets. The basic dif¬ 
ference is that magnetic systems are neutral Bose gases El, 
while the charged Cooper pairs are coupled to the electromag¬ 
netic gauge field. In other words, the orbital coupling to an 
external field that stabilizes the Abrikosov vortex crystal in 
type-II superconductors is basically absent in magnets. 

Topological spin structures must then be stabilized by other 
means. A key aspect of magnetic systems is that compet¬ 
ing interactions are ubiquitous. A common outcome of this 
competition is a magnetic susceptibility that is maximized by 
several low-symmetry wave vectors Q connected by point- 
group transformations of the underlying material. Topological 
spin structures can emerge when the effective interaction be¬ 
tween the different Q modes favors a multi-g ordering. This 
is the case of the 520 materials, in which the Dzyaloshinskii- 
Moriya inoiini interaction D that arises from their noncen- 
trosymmetric nature shifts the susceptibility maximum from 
Q - Q favored by the ferromagnetic exchange 7 to a finite 
vector Igl DIJ that can have different orientations due to 
the cubic symmetry of the 520 structure. Thermal fluctua¬ 
tions then play an important role for stabilizing the 6-Q struc¬ 
ture that leads to the hexagonal skyrmion crystals in bulk ver¬ 
sions of the B2Q materials 0. In contrast, the phase is already 
stable at the mean field level in 2D thin films El- In addi¬ 
tion to chiral magnets, skyrmion crystals lfT2l . soliton crys¬ 
tals lfT3l [T4ll . and Z2 vortex crystals ifTSll have been theoret¬ 
ically predicted in other classical spin systems. All of these 
examples correspond to 2D crystals of topological structures, 
i.e., they are not modulated along the third dimension. 

More recently, two of us proposed the realization of mag¬ 
netic vortex crystals in a quantum spin system of weakly cou¬ 
pled triangular layers near a magnetic field-induced quantum 
critical point (QCP) ifThll . The basic idea is to use geometric 


frustration as the source of competing interactions and quan¬ 
tum fluctuations to stabilize the multi-g vortex crystal states. 
This study focuses on a case with six degenerate coplanar Q 
vectors that are connected by the Ce symmetry transforma¬ 
tions of the underlying lattice. Consequently, as in the pre¬ 
vious examples, the resulting vortex crystal is not modulated 
along the third direction. 

In this Letter, we demonstrate that a similar mechanism can 
also stabilize exotic 3D crystals of vortex lines. Unlike the 
case of the 2D vortex crystals, we are unaware of alternative 
realizations of 3D vortex crystals. As we explained above, the 
observation of magnetic skyrmion lattices unveiled the rele¬ 
vance of multi-g orderings that produce 2D crystals of topo¬ 
logical structures. However, much less effort has been devoted 
to the 3D crystals that can also arise from multi-g orderings. 
The recent real-space observation of a skyrmion-antiskyrmion 
cubic lattice in MnGe 03 confirms the physical relevance 
of these 3D structures. The key to realize 3D crystals of 
topological objects is to find regions of stability of multiple 
noncoplanar-Q orderings. Consequently, we study the body- 
centered (bcc) and face-centered cubic (fee) lattices that com¬ 
monly occur in nature (typical examples are the transition- 
metal oxides and fluorides ifTSl . solid ^He ifTOll . 3D Wigner 
crystals ll20l . and the alkali-metal fulleride CssCeo EIIES)- 
By extending the exchange interactions up to third nearest 
neighbors, we produce a single-magnon dispersion with mul¬ 
tiple degenerate minima at noncoplanar Q vectors connected 
by the cubic point group. We compute the optimal single¬ 
particle state for condensing the magnons and find that several 
multi-g states corresponding to different vortex crystals span 
sizable regions of the phase diagrams with isotropic exchange. 
These phases are further stabilized by symmetric exchange 
anisotropy that arises from, e.g., dipole-dipole interactions or 
spin-orbit coupling. The resulting spin textures consist of ex¬ 
otic 3D patterns of vortex strings. 

We consider a spin-5 Heisenberg model on bcc and fee lat¬ 
tices coupled to a magnetic field; 
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FIG. 1. (color online), (a)-(f) Schematic momentum-space represen¬ 
tations of the multi-g condensates near the field-induced QCP for the 
case of six degenerate minima. The arrows representing the phases 
On of the g„ component of the order parameter [see Eq. ([^] are only 
shown for states in which their relative values are fixed by the inter¬ 
actions or anisotropy. The gray (light) color indicates no correlation 
among the different phases 0„. 


where 7^ are the Heisenberg interactions up to 3rd nearest 
neighbor {Jx,J 2 ,Ji]. In this study, we focus on the external 
field H applied along the high symmetry [111] direction. The 
spin -2 operators can be represented by hard-core bosons ^23\ : 
Sf - hi, S~ - b], S^. - 112 - b\bi, where the z axis is along 
the magnetic field direction. The Hamiltonian is thus trans¬ 
formed into a model for an interacting Bose gas: 

^ Yj^U + yq)bl^qbl,^^bk'bk, (2) 

h k,k\q 

where Wfc is the single-boson (magnon) dispersion, /t = //sat - 
H is the chemical potential, // ^ oo is the on-site hard-core 
potential, and Vq (Fourier transform of Jij) is the density- 
density interaction arising from the Ising component of the 
spin exchange Il24l . 

The relative strengths of /i, J 2 , and Jj, determine the num¬ 
ber of degenerate minima in the single-magnon dispersion Wfc. 
Phases with six and eight minima exist in both the bcc and 
fee lattices. A phase with twelve minima also exists in bcc 
lattice El. For concreteness, we will focus on the region 
with six degenerate minima, whose positions are denoted by 
±Qn - ±Ge,j, where n = 1,2,3. 

The single-magnon dispersion becomes gapless at // = //sat 
which signals the phase transition into a Bose-Einstein con¬ 
densate IfThl I25l - l3^ . In the vicinity of this transition \H\ < 
|//satl, the boson density is vanishingly small, and we can use 
Beliaev’s dilute boson approach ll^ to compute the effective 
boson-boson interactions in the long-wavelength limit. Be¬ 
cause this is a controlled expansion in the small lattice gas pa¬ 
rameter (ratio between the scattering length and the average 
interparticle distance), the result is asymptotically exact in the 
dilute limit. The next step is to condense the bosons in the 
most general single-particle state, i.e., to replace the bosonic 
operators for each wave vector Q by six complex amplitudes: 
ib±Qjl 'JN - ^Jp±Q^ exp J. The total energy is the sum 
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FIG. 2. (color online). Phase diagrams of the Heisenberg model 
Q under the nearly saturated magnetic field, where Wfe has six de¬ 
generate minima in the colored phases, (a) bcc lattice, /, > 0, no 
anisotropy. “NAF” denotes the case with Q = (2n, 2n, 2n). (b) bcc 
lattice, /i > 0, anisotropy dominating region, (c) bcc lattice, /i < 0, 
without anisotropy. “Ferro” denotes Q at (0,0,0), and “PS” denotes 
regions where we have phase separation or bound states, (d) bcc lat¬ 
tice, /] < 0, anisotropy dominating region, (e) fee lattice, 7, > 0, no 
anisotropy, (f) fee lattice, anisotropy dominating region. 


of the low-energy terms allowed by translation symmetry, i.e., 
density-density interactions between bosons in the same (F 1 ) 
and different (F 2 ,F 3 ) minima, as well as a F 4 vertex that scat¬ 
ters bosons between two pairs of opposite minima, 

^ ^ y ^ PqP-q„ +^ 3 X 1 

n,(T=± n n<m 

(T\ ,cr2=± 

^'^^^YjJpqP-qPq.P-Q,,, - PP, ( 3 ) 

n<m ~ 

where p — Yjn {pQ„ + P-e„) the total boson density and 
*hn = ‘pQ„ + 4>-Q„- The interaction vertices Fi,... ,F 4 are ob¬ 
tained by summing over the ladder diagrams at zero total fre¬ 
quency 12^ . 

The zero-temperature phase diagram is determined by min¬ 
imizing the total energy E given in Eq. QUm). Depending 
on the relative strengths of exchange interactions, one of the 
six possibilities in Eig.[T]is realized. Out of these, three con¬ 
densates in particular realize vortex crystals: 3-Q I, 4-Q I, 
and 6-Q 11 (Fig. [^. Another reason for considering these 
states is that the latter two are further stabilized by symmet- 
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lie exchange anisotropy originated from spin-orbit coupling 
or dipole-dipole interactions. Close to the saturation field //sat, 
this exchange anisotropy yields the interaction term ll24l : 

Eaoc JA yjPQ„P-Q„ cos(<t)„ + 2nnl3 - nil). (4) 

n 

Although Ja is typically small, Ea ~ \p\ is linear in the boson 
density. Consequently, in the dilute limit {p <Si \Ja/Ji\ 1), 
it always dominates over the exchange interaction in Eq. 0. 
The calculations show that the 2-Q I, 4-Q I, and 6-Q II con¬ 
densates are the three lowest energy states, degenerate to lin¬ 
ear order in the density p. The degeneracy is lifted by further 
considering the second-order density-density interactions in 
Eq. @, stabilizing the 6-Q II state over a wide range of pa¬ 
rameters on both bcc and fee lattices (see Fig.|^. On the other 
hand, sufficiently far away from the QCP (but still in the low 
density regime, \JaIJi\ p ^ Ea is negligible andEq. 0 
alone determines the ground state configuration. The result¬ 
ing phase diagrams for negligible (dominant) anisotropy are 
shown in the left (right) column of Fig. respectively. 

We now focus on the 3-Q I, 4-Q I, and 6-Q II states that 
realize vortex crystals whose spin structures on [111] layers 
are illustrated in Fig. We find that the vortex and antivor¬ 
tex cores form regular lattices on every layer. Vortices and 
antivortices correspond to different signs a: = +1 of the vec¬ 
tor spin chirality S', x Sj when we circulate (i j) around 
the vortex core. Explicitly, in the vicinity of the (anti)vortex 
core arctan(5]'/5p = arg{/>|) ~ K<p -i- y 6{ip), where y is the 
helicity |[35l and 6{(f) is a 27r-periodic function of the polar 
angle ip around the (anti)vortex core in the [111] plane such 
that d(p6{p) - 0 and |i5(i/j)| <si 27 t. For a given chirality, the 
(anti)vortices can have different relative helicities, as we show 
in Fig. the 3-Q I state includes three types of vortices with 
helicities that differ by 27r/3 in each [111] plane (the same is 
also true for the antivortices), see Fig.[^; similarly, two types 
of (anti)vortices appear in the 4-Q I state with helicities that 
differ by n [Fig. |^]; the 6-Q II state contains three types of 
vortices with helicities that differ by 27r/3 and two types of 
antivortices with helicities that differ by n [Fig.j^]. 

The vortex cores are strings that extend along the third di¬ 
mension. These strings form different patterns for each multi- 
Q condensate. The vortex and antivortex strings form parallel 
straight lines along the [111] direction in the 3-^1 and 4-Q 
I states. The same is true for the antivortices of the 6-Q II 
states. However, the vortex strings form a more exotic pat¬ 
tern in the 6-Q II state. As is shown in Fig. the vortex 
strings cross each other and the helicities of the crossing vor¬ 
tices and antivortices are shifted by n. This unusual behavior 
arises from the fact that the six Q vectors are noncoplanar. In 
contrast, when the condensate Q vectors are on the same plane 
in the reciprocal space, as is the case with the 3-Q I and 4-Q 
I states, and those considered in Ref. ea, the vortex strings 
are straight lines along the high-symmetry axis. The helicity 
of each (anti)vortex increases linearly in the layer index for 
the 3-Q I state. In contrast, the helicity of each (anti)vortex 
is shifted by tt between consecutive layers of the 4-Q I state. 



FIG. 3. (color online). Spin structures in the vortex crystal phases 
on [111] layers (H || [111]); only the xy components are plotted. 
The color intensity denotes the boson density p (the transverse spin 
components are ~ yfp) and the bright spots denote fully polarized 
spins. The circles with (without) crosses denote vortex (antivortex) 
cores. Different colors of the circles denote different helicities y (see 
text), (a) The 3-Q I state on the bcc lattice for |g„| <K 1. (b) The 4-Q 
I state on the bcc lattice for |g„| ;£ In. (c) The 6-Q II state on the bcc 
lattice for |g„| <K 1. 

This alternation arises from the fact that |g„| < 2ji and each 
[111] layer belongs to a sublattice of the bcc lattice. In the 6- 
Q II state, the helicity of each (anti)vortex is constant except 
for crossing points where it is shifted by n. 

While in this Fetter we only focus on the regions where 
(i)k has six degenerate minima, there are other regions in the 
phase diagrams where vortex crystals should also arise. For 
example, u)k with eight minima is realized in a region next to 
the 6-minimum region on both the bcc and fee lattices, and 
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FIG. 4. (color), (a)-(d) The sequence of layers along the [111] di¬ 
rection for the 6-Q II state on the bcc lattice in the case \QJ <K 1. 
Several intervening layers between (a) and (b), etc., are omitted for 
simplicity. The arrows indicate how the vortex cores (circles with 
crosses) move from one layer to the next layer above. The lines of 
the antivortex cores (circles without crosses) are parallel to the [111] 
direction. The region enclosed by a square follows one of the an¬ 
tivortex core. The n helicity shift occurs below the layer (c). (e) 
Three-dimensional picture of the vortex string structure. The strings 
corresponding to antivortices in (b)-(d) are not shown. 

a 12-tninimutn case also occurs on the bcc lattice El- A 
similar calculation can be applied to these cases in order to 
obtain the stable multi-g orderings. 

We note that when some of the effective interactions are 
attractive (typically the case with ferromagnetic exchanges), 
the system may undergo a hrst order phase transition at the 
saturation held il3^ [^ . This implies that it is unstable to¬ 
wards phase separation if one hxes the particle number (i.e., 
the z component of magnetization) in the canonical ensemble. 
An alternative scenario is a continuous transition associated 
with the condensation of multimagnon bound states We 
have verihed that none of these cases takes place for antifer¬ 


romagnetic nearest neighbor interactions (Ji > 0) in both the 
bcc and fee lattices. For ferromagnetic nearest neighbor in¬ 
teractions (7i < 0), on the other hand, the phase separation 
occurs in a large region of the phase diagram, as indicated in 
Figs.l^c) and|gd). 

Although crossings of vortex lines have been observed in 
superconducting vortex glasses and liquids ifJTl [38], we are 
unaware of the existence of 3D vortex crystals like the one 
shown in Fig. These vortex crystals can be detected with 
neutron diffraction in single-domain samples. Materials with 
more than two degenerate minima in the magnon dispersion of 
the fully saturated state can be identihed directly by measur¬ 
ing the inelastic neutron scattering spectrum at |//| > |//satL or 
indirectly, by extracting the exchange constants from the zero- 
field inelastic neutron scattering spectrum. Nuclear magnetic 
resonance (NMR) also allows us to distinguish among differ¬ 
ent multi-g orderings because the NMR line shape is in gen¬ 
eral qualitatively different for single-, double-, and three-g 
orderings. 

Our results indicate that these materials are strong can¬ 
didates to exhibit magnetic vortex crystals just below their 
saturation field. While here we have considered the partic¬ 
ular cases of bcc and fee lattices as examples, the general 
principle can be directly extended to other highly frustrated 
structures, such as hep and pyrochlore lattices, which are 
also common in nature. Based on our calculations, it is ex¬ 
pected that exchange anisotropy (due to e.g., dipolar interac¬ 
tions) will select a double-g magnetic ordering or a magnetic 
vortex crystal. The selection mechanism between these two 
competing phases is provided by the effective interaction be¬ 
tween magnons, which ultimately depends on the details of 
the exchange couplings. There are several candidate materials 
that comprise highly frustrated 3D lattices of rare-earth mag¬ 
netic ions. Because the exchange anisotropy is expected to 
be stronger in these ions due to a large spin-orbit interaction, 
the vortex crystal phase could extend over a wider window of 
magnetic held values. 
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TRANSFORMATION FROM THE SPIN TO THE BOSONIC LANGUAGE 


—Supplemental Material— 


Transformation from the spin to the bosonic language 


The spin- 5 Heisenberg model is defined on the bcc and fee 
lattiees, shown in Fig. SI By using the Matsubara-Matsuda 
transformation introdueed in the main text, the Hamiltonian is 
transformed into the bosonie language, up to a eonstant term: 




Jl 


(ij) 




((<7» 

U 


((('7»> 


-I- 7i ^ riiiij + nitij + Ji'^ nirij 

(ij) ((ij)) (((ij))) 

(Z\J\ + Z 2 J 2 + Z 3 J 3 




(SI) 


where U is the on-site hard-eore repulsion, whieh is sent to 
infinity in the ealeulation, and Zi,Z2,Z2 are the eoordination 
numbers of the 1 st, 2 nd and 3 rd nearest neighbors. 


(a) bcc Jl > 0 
1 


0.5 


(b) bcc Jl < 0 


G-min 

12-min 

(27r, 271,27r 

) 

iij,-'— 


-1 

(c) fee Jl > 0 
1 


6-inin 

12-min 

(0,0,0 

) 



(d) fee Jl < 0 



6-min 



( 0 , 0 , 0 ) 


-1 - 0.5 0 0.5 1 -1 - 0.5 0 0.5 1 

72/17i| 72/|7i 


FIG. S2. The single-magnon phase diagrams, with each region de¬ 
noted either by the number of minima (when there are multiple min¬ 
ima at incommensurate Q-vectors), or denoted by the positions of the 
Q-vectors (when Q is commensurate). (a)(b) bcc lattice. (c)(d) fee 
lattice. 




FIG. Si. The Heisenberg interactions Ji,J 2 ,Jt, are defined on the 
1st, 2nd and 3rd nearest neighbors, (a) bcc lattice, (b) fee lattice. 

By Fourier transformation b\ = the 

Hamiltonian is written down in fc-space: 

H^Yj^e{k)-e{0) + H]blbk 

k 

^ fc.fc'.q 


where 



m m Tj3 

here r,j denote the positions of the neighboring sites. And 


For fee lattiee: 


e(k) - 2Ji(^ eos — eos — -i- eos — eos — 
kx \ / , 

+ eos — eos — j + J2{^ eos kx + eos ky+ 

\ / k^ 

eos + 4J3(^ eos kx eos — eos —h 


2 

kx j kx ky >. 

eos kx eos — eos-h eos k, eos — eos — 

■ 2 2 2 2 / 


(S6) 


We define the minimum value of e(k) to be e^in, in this 
way Wfc s e{k) - en,i„ has minimum value equals to zero. The 
Hamiltonian is rewritten as: 

-E)blbk + ^ ^(t/ + Vq)bl^y^,^^bk'bk (S7) 

k k,k',q 

where the ehemieal potential: 

fi = [e(0) - e,nin] -H = H,,r-H (S8) 


Vq = 2 e{q) (S 4 ) 

To be explieit. For bee lattiee: 

kX by k^ / 

e(k) - 4 Ji eos — eos — eos — + J2y eos kx 

+ eos ky + eos k^ ) + eos kx eos ky 

+ eos ky eos k^, + eos eos k^ (S 5 ) 


Beeause of the frustration, the single magnon dispersion ojf, 
ean have multiple degenerate minima at different g-veetors. 
In Fig.|S 2 [ we eompute the number of minima in Uk, for both 
bee and fee lattiees. 

For eonereteness, we foeus on the regions with 6 degenerate 
minima, whose positions are denoted by +g„ = ±ge„, where 
n- 1 , 2 , 3 . The value of g is given by eos ^ =-Ji/(J2-1-4J3) 
for the bee lattiee and eos ^ = -(Ji + 2 J 3 )/(J 2 -1- 4J3) for the 
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EFFECT OF SYMMETRIC EXCHANGE ANISOTROPY 


fee lattiee. Correspondingly, the saturation field values are: 
2/2 

+ 272 + 873 (S9a) 

72 + 473 

+ 47, + 272 + 1673 (S9b) 

72 + 473 

Calculation of Effective Interactions 


simplicity, and (F) = f ^^Fg/. We also assume that Vg is 
centro-symmetric, ie 

J d^gV(g)^0 (S13) 

By substituting the ansatz into the Bethe-Salpeter equation 
and taking the hard-core limit, we get the following form of 
linear equations: 


The effective interactions in the dilute limit for hard-core 
bosons are calculated by the Bethe-Salpeter equation, which is 
equivalent to summing over all the ladder diagrams (Fig. S31. 


rg(k,k')= U + Vg 


f 


CPg' Tg,{k,k'){U +Vg^g.) 
Fez Wfc+q' -H Uk'-g' 


(SIO) 


+ 6gy)A, + tKD = 1 

V 

where the integrals are defined as: 


(S14a) 

(S14b) 


where Fez is the volume of the 1st BZ. 



To = ^ 

r cPq 1 

I Fez Wfc+g + a>k'-q 

t'' - 

r (pq 


I Fez Wfc+g + a>k'-q 


r Pq 

^ j 

1 Fez ^k+q + 


(S15a) 

(S15b) 

(S15c) 


FIG. S3. Ladder diagrams. 


When the magnetic field H is close to the saturation value 
i/sat, the system is unstable towards BEC at the dispersion 
minima. In this case we can take the long wave length limit 
k —> ±Qi, and calculate the corresponding vertex functions 
(schematically shown in Fig.|S4|i: 


r, = To{Q„,Qn) 

Ti - ro(g„, -Q„) + F_2e„(6„, -Qn) 

r3 = ro(g„, Q,„) + Tq_^-q^(q„, g„,) 

^4 = Tq (g„, -g„) + F_q _g (g„, -g„) 



Denote: 


Bgy = TJ^ViVy) + 6gy (S16a) 

Cg = V{rg){T\r (S16b) 


The above equations are now organized into a matrix form: 


'fill ■ 

• fiu tJ ' 


■ 


f ^ 

fi.l ■ 
Cl ■ 

■ t \ 

■ Cj To . 


A. 

Ur) J 


1 

. 1 j 


(S17) 


By solving the linear equations Eq. ( |S17[ l, we obtain all the 
unknown coefficients in the ansatz Eq. ( S12[ ). Then we can 
substitute the values of F, ,..., F 4 into the expression of effec¬ 
tive energy, and determine which multi-g state will be stabi¬ 
lized. 


Effect of symmetric exchange anisotropy 


FIG. S4. Schematic plot of the vertex functions in the long wave 
length limit. 

To solve the Bethe-Salpeter equation, we start from the fol¬ 
lowing ansatz: 

Tg = <F) + ^A,F(r,y'J'-'' (S 12 ) 

n 

where denotes the positions of the 1st, 2nd, and 3rd neigh¬ 
boring sites. Tthe k, k' index in Fg(fc, k') are omitted for 


We consider short-range symmetric exchange anisotropy 
(cutoff at 2 nd nearest neighbor): 

HA^J]-3(Srrij)(Sj-rij) (SI 8) 

(U) 

such terms can arise directly from dipole-dipole interactions, 
or perburbatively from spin-orbit coupling lfTTl . 

Similar to the treatment of the Heisenberg exchange inter¬ 
actions, we choose the quantization axis along [ 111 ] direction. 
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EFFECT OF SYMMETRIC EXCHANGE ANISOTROPY 


and represent the spin -5 operators with hard-core bosons. In 
the long-wavelength limit, for both bcc and fee lattices; 


Ha oc 


V3 


V3 






2' G 2 -Qi 




-i- h.c 


(S19) 


Then we condense the bosons by {b±Q )l y/N - 
■\jp±Q^ exp (i(p±Q^^, which gives the energy correction of sym¬ 
metric exchange anisotropy: 

Ea^ Ja'^ ylPQ„P-Q„ cos(0„ -H 2nnl3 - njl). (S20) 


where O,, = (pQ^ + <f>-Q„. 
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